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ABSTRACT 

Active galactic nuclei are clearly heating gas in 'cooling flows'. The effectiveness and 
spatial distribution of the heating are controversial. We use three-dimensional simula- 
tions on adaptive grids to study the impact on a cooling flow of weak, subrelativistic 
jets. The simulations show cavities and vortex rings as in the observations. The cavi- 
ties are fast-expanding dynamical objects rather than buoyant bubbles as previously 
modelled, but shocks still remain extremely hard to detect with X-rays. At late times 
the cavities turn into overdensities that strongly excite the cluster's g-modes. These 
modes damp on a long timescale. Radial mixing is shown to be an important phe- 
nomenon, but the jets weaken the metallicity gradient only very near the centre. The 
central entropy density is modestly increased by the jets. We use a novel algorithm to 
impose the jets on the simulations. 
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1 INTRODUCTION 

The arrival of data from the Chandra and XMM-Newton 
X-ray satellites has initiated a major change in our un- 
derstanding of cooling flow systems. In particular, the new 
data have convinced many astronomers that these sys- 
tems are powerfully influenced by AGN. Many aspects 
of this interaction are controversial, however, because the 
dynamics of the interaction between high-entropy plasma 
ejected by the AGN and the surrounding thermal plasma 
is complex. In particular, the interaction is inevitably time- 
dependent, chaotic and geometrically irregular. Numerical 
simulations are capable of giving us insights into the in- 
teraction and current thinking has been strong ly influenced 
by the simulations of ICh urazove^aUj^OfJlJ^IOirnise^alJ 
l200j) . lBriiggen fc Kaisd feOOlLlBriiggen fc Kaiserl |2002>) . 
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Briiggen et al. I feOCffl. iRevnolds et alj fcOOll |2002I> . and 
Basson fc AlexandeJ i200ot) . In this paper we attempt to 
clarify the issues and to explore new territory with numerical 
simulations that are novel in several respects. Most crucially 
we present the first three-dimensional simulations that use 
an adaptive grid to simulate the entire cooling flow with sig- 
nificant spatial resolution at the centre. We draw attention 
to the possible importance of momentum input from the 
AGN for the dynamics of the cluster gas, and explore a part 
of the parameter space that is opened up once momentum 
input is considered. 

In Section 2 we motivate our simulations by reviewing 
the history of cooling flows and the reasons why the new data 
have had a profound impact. We explain why simulations 
are essential for developing a sound understanding of the 



connection between a cooling flow and its embedded AGN. 
In Section 3 we motivate our phenomenological model of 
how an AGN impacts on its immediate surroundings. In 
Section 4 we introduce a new way to simulate jets that are 
embedded in a large-scale system. Section 5 describes a pair 
of simulations, and Section 6 compares these simulations 
with previously published ones. In Section 7 we discuss the 
implications of our work and what should be done next. 



2 COOLING-FLOW HISTORY 

The early X-ray satellites discovered that the deep gravita- 
tional potential wells of clusters of galaxies and luminous el- 
liptical galaxi es confine substantia l masses of gas at the virial 
temperature iF orman et al.||l972|: Bahcall & Bahcall|JlJ)7f& 
Mitchell^^dJ ^7^: [Se^temitsos e^lj ^)77]: [Formaiie t alJ 
1978: lAbramopoulos fc Kulll983ilJones fc Formanlll984l) . It 
was recognized that towards the centres of th ese system s 
the cooling time is shorter tha n the Hubble time 
Motivated by this observation. ICowie fc Binnevl 1(1977) mod- 
elled these systems under the assumption that they had 
achieved a steady state, in which gas seeped slowly in- 
wards as it radiated energy. The early models did not ex- 
plain what became of gas once it had cooled to below the 
virial temperature and out of the X-ray band, beyond argu- 
ing that isobaric cooling is thermally unstable, so at suffi- 
ciently _smdl_j^adfr_clou(fe to 
form dCowie fc Binnevlll977D . iFabian fc Nulsenl (119771) ar- 
gued that these cold clouds gave rise to the Hcv filaments 
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that were observe d in many cooling flows jLvndsl Il97(il 
ICowie et al."lll983l) . 

As the spatial resolution of the X-ray spectroscopy in- 
creased, it became clear that if the steady-state conjecture 
was valid, mass must be dropping out of the flow into cool 
clouds throughout the region of radius r < r coo i, where r coo i 
is the ra dius at which the cooling time equals the Hub- 
ble time. iNulsenl 1^)86) introduced a theory of a multi- 
component intracluster medium (ICM) which allowed for 
such distributed 'mass dropout', and showed that the ob- 
served X-ray brightness profiles of cooling flows, which were 
less steep than the original steady-state models predicted, 
could be explained if the distribution in temperature of the 
gas at a given radius were suitably chosen. 

A puzzling feature of models with distributed mass 
dropout was the failure of strenuous efforts to observe ob- 
jects - ionized or neutral clouds, or young stars - that 
formed from the 10 to lOOOMp yr _1 of dropout that the 
models predicted ]Fabianlll994l). Ionized gas was detected 
early on jLvnddmfflT l Cowie et al. I Il983|) . and molecular 
gas has been detected more recently ]Doiiahue et alJl2000t 
lEdgdl200ll: lEdee et alJl2002fl . but the mass of such matter 
falls short of predictions by a factor of at least ten, and it is 
more centrally concentrated than predicted. As the quality 
of the available spectra improved, the models' predictions for 
the X-ray spectra also came into conflict with the data. In 
particular, the X-ray luminosity be low 1 keV, tended to be 
significantly smaller than predicted JStewart et alJl984bllal) . 

2.1 Cracks in the edifice 

Soon after the tenth birthday of cooling-flow theory its 
intell ectual foundatio ns were shot away by two papers. 
iMalaeoli et alJ <ll987t) showed that when gas at a range of 
specific entropies is confined by a potential well with specific 
entropy increasing outwards, the gas is not thermally unsta- 
ble: a region in which the temperature is lower than ambient 
is under-buoyant and will sink and then oscillate around the 
radius at which the ambient medium has the same specific 
entropy. In pr inciple, cooling causes these oscillati ons to be 
over-damped fealbus fc Sokejll989t lTribblelll989ll . but the 
growth rate is low and in a real system turbulent damping is 
likely to suppress the over-stability. In a s econd lethal attack 
on th e steady-state cooling-flow model, iMurrav fc Balbusl 
( 1992) showed that the differential equations that govern 
trapped cooling gas do not tend to a steady state. Rather, 
as the system evolves from a generic initial condition, the 
importance of partial derivatives with respect to time in- 
creases secularly because the s ize of dx/dt is no t of or der x 
over the age of the system, as lCowie fc Binnevl 1119771) and 
subsequent papers had implicitly argued, but ~ x/(t c — t), 
where t c is the time of a future cooling catastrophe, when 
the local density will diverge. 

In view of these theoretical developments, a handful 
of theorists argued that a cooling flow is an unsteady re- 
sponse to heating by a galactic nucleus, the nuclear activ- 
ity being episodically stimulated b y the development of a 
cooling catastrophe in its environs jTabor fc Binnevlll993t 
Binnev fc TaboJl995llBinnevll99^ICiotti fc Ostrikeiil997i 
Binnevl Il999t ICiotti fc OstrikeJ l200lD . Unfortunately, the 
steady-state cooling-flow picture had by this stage devel- 
oped such momentum that it could not be derailed by mere 



theoretical considerations. A variety of mechanisms was in- 
voked to explain the conflict between the basic model and 
the new data. These included strong internal absorption 
]Allen fc Fabian lll997f). magnetic l ocking of over-dense re- 
gions jTribbld Il989t balbusl Il99lt) abundance anomalies 
( Fabia n et al . 200 1; Mo rris fc Fabianll2003ft. and conduction 
of heat from large to small radi i dBertschinger fc Meiksinl 
119861: iNaravan fc Medvedevll200ll) . 

2.2 Decisive new data 

The arrival of data from Chandra and XMM-Newton has 
finally turned the tide of opinion away from the dis- 
tributed mass-dropout hypothesis. In particular, the new 
data show that internal absorption is not viable as an ex- 
planation of the paucity of flux below 1 keV because it sup- 
presses very soft phot ons that are observed in abundance 
ilBohringer et all2002l) . Moreover, spatially resolved spectra 
are incompatible with the predicted multiphase ICM: there 
is no evidence fo r plasma at more than one temperature a t 
a given location dTamura et aljl200ll : iPeterson et al. ll2002l) . 
Within the cooling radius the temperature of the ICM is 
found to decline with decreasing radius, but has a floor value 
that varies from system to system. The floor almost always 
lies above 1 keV, so the paucity of flux below 1 keV simply 
reflects an absence of very cool gas, in clear conflict with the 
steady-state model. [Kaiser fc Binnevl |2003i) have shown, by 
contrast, that if cooling flows are episodically heated, very 
cool gas will rarely be seen, because in any given system it 
appears only fleetingly and in small quantity. 

In addition to completing the falsification of the mass- 
dropout model, the new data have provided morphologi- 
cal evidence for heating by AGN. It has long been rec- 
ognized that cooling f lows generally have e mbedde d non - 
thermal radio sources jFabbiano. Gioia fc Trinchierilll989T) . 
but it used to be argued that there was no convincing ev- 
idence that th ese sources had a significant impact o n their 
cooling flows 1 Allen et al.l l2001i iFabian et alJ 1200 IT), even 
thoug h the beautiful radio map of Virgo A bv lOwen et alJ 
(2000) showed that in Virgo AGN-powered ultrarelativistic 
electrons densely permeate the region r < r coo \/2. Imag- 
ing data for several systems now show 'cavities' in the 
X-ray emitting plasma that coincide to s ome degree with 
peak s in the non-thermal radio emission feohringer et alJ 
[ 19931: IChurazov et al] l200d : iMcNamara et al] 12001 l200ll : 
blanton et all200lL The interpretation that the cavities are 
'bubbles' inflated by r adio sources in the the r mal plasma, 
has gained currency (lClTurazov_el_alJ |200jJ; |Quflis_et_al] 
120011: iBriiggen fc Kaiserl 1200 ll 120021: iBriiggen et al. 1120021) . 
Consequently, it is now widely believed that cooling-flow gas 
is significantly affected by AGN. 



3 KEY QUESTIONS 

Many fundamental questions regarding the nature of 
AGN/cooling-flow interaction remain to be answered. For 
example: 

(i) Are the observed cavities inflated by ultrarelativis- 
tic jets or by sub-relativistic bipolar flows, or by a com- 
bination of the two? We shall see that the answer to this 
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question has a strong bearing on the dynamics of the cav- 
ities: an ultrarelativistic jet carries very little momentum, 
so the cavities it inflates simp l y rise buoyantly, as in th e 
simulations of [ C hurazo v et alJ i200ll) . iQuilis et al] feOOll) . 
iBriiggen fe Kais*eiri2002ll and others: a subrelativistic wind 
imparts momentum as well as energy to the cavity, which is 
thus driven up through the ambient gas, with implications 
for the stability and longevity of the cavity. 

(ii) How is the energy that we know lies within the ob- 
served cavities, transferred to the thermal X-ray emitting 
gas? 

(iii) Do the observed cavities in cooling flows such as 
those in Hydra and Perseus heat the systems fast enough to 
offset radiative cooling? While it is now clear that cooling 
flows are heated by AGN and that locally partial derivatives 
with respect to time are significant terms in the governing 
equations, it is not clear whether the azimuthally averaged 
profile s of flows are strongly time dependent. At one ex- 
treme iKaiser fc Binnevl <l2003f) present a picture in which 
heating is unimportant for several hundred megayears be- 
tween catastrophic AGN outbursts, and the radial density 
profile evolves significantly between outbursts. At the other 
extreme there is the possibility that the cores are in a tur- 
bulent steady-state in whi ch azimuthally-average d heating 
and cooling are in balance llTabor fc Binnevlll993l) . 

(iv) On what timescale does a cooling-flow settle to ap- 
proximate hydrostatic equilibrium after a nuclear outburst? 

(v) How visible in X-rays will be non-equilibrium struc- 
tures, such as cavities and shocks, that the AGN generates? 

(vi) What is the radial distribution of injected energy af- 
ter approximate hydrostatic equilibrium has been restored? 
This distribution effectively determines the predicted X-ray 
brightness profile and the time elapse between outbursts. 

(vii) What impact does AGN-induced radial mixing have 
on metallicity gradients? Is this impact compatible with 
models of metal enrichment and observed metallicity dis- 
tributions? 

(viii) We must expect the outflow from the AGN to fluc- 
tuate on very short timescales. What observable phenom- 
ena will be generated by these fluctuations? A priori we ex- 
pect the shortest timescale fluctuations to be evident at the 
smallest radii. What phenomena are expected at different 
radii, and what are the associated timescales? 

(ix) Jets from accreting objects ar e known to pre cess - 
most famously in the case of SS 433 <Milgromlll97sh . How 
would the inflation of cavities be affe cted by jet p recession 
or by a wide opening angle of the jet (Sokcr 2003)? 

(x) Gas trapped in a potential well is likely to have non- 
zero angular momentum. As the gas cools, the dynamical 
importance of angular momentum is liable to increase - in 
proto-spiral galaxies it evidently becomes dominant. Radial 
stirring of gas by an AGN is liable to move angular momen- 
tum outwards and make it dynamically less important. How 
rapid is such angular-momentum transport? What impact 
does it have on the evid ent failure of giant ellip tical galaxies 
to form disks recently femsellem et al.|l2002f) ? and on the 
locati on of the molecular gas that has been found in cool ing 
flows JPonahue et alJl2000l : lEdgeil200ll : TEdee et al.ll2002l) ? 

Since the interaction between an AGN and a 
cooling flow is inherently non-spherical and unsteady, 
simulations have a vital role to play in answering 



these questions. Several sets of simulations of AGN / 
cooling-flow interaction have appeared in the literature 
jjCjmrazov ej i _alj|2j)0jyQvrUig_et aj]|200iyPmij ^^^ 
| 200lll2002l:lRevnolds et alJ200lll2002l:lBruggen et al. 120021 : 

basson fc Alexander! 120031) . For reasons of computational 
cost, most published simulations assume spherical or ax- 
isymmetric symmetry or use a rather coarse spatial resolu- 
tion. We present simulations that are fully three-dimensional 
and achieve competitive spatial resolution in a large box by 
exploiting adaptive grids. 

3.1 Heating mechanism 

Published simulations adopt a variety of approaches to the 
way in which the AGN t ransfe rs energ y to its environs. 

ICiotti fc Ostrikerl (1 1 997t 120011) consider inverse- 
Compton heating of the gas. This is an inefficient process 
since a photon of energy E transfers to the gas a fraction 
~ E/m. e c 2 of its energy, so even if the Thompson-Compton 
optical depth to the AGN is significant, radiation carries 
away most of the AGN's output unless the spectrum of the 
radiation is hard. Moreover, for observed cooling flows the 
Thompson-Compton optical depth from infinity to radius 
r is small down to the smallest radii (r > 100 pc) that can 
be resolved in state-of-the-art simulations. Consequently, 
any inverse-Compton heating of the gas is likely to be 
concentrated at unresolved radii, and to impact on the 
simulated region by forcing an inflow across the simulation's 
inner boundary. In this respect the inverse-Compton model 
is similar to the model on which we shall concentrate. 
The essential difference is that we shall assume a strongly 
bipolar outflow, whereas inverse-Compton heating might 
produce a nearly spherical flow at the inner boundary. 

Our simulations assume that a cooling flow is primar- 
ily heated by a bipolar flow. This assumption is based on 
several considerations. First it is now generally accepted 
that AGN are powered by accretion onto a massive black 
hole. Bipolar outflows are obs erved around most accr eting 
objects, be they forming stars jArce fc Goodm^[200j ), ac- 
crcting black holes that power AGN JPounds et al. 1120031) 
and micr o-quasars jMirabel fc Rodriguezlll999l) . or forming 
galaxies iPettini 61^012*0021) . Hence, it is natural to expect 
a bipolar outflow to emerge from an AGN that is heating a 
cooling flow. 

A second reason to believe that cooling flows are 
heated by bipolar flows follows from the observation 
iFabian fc Canizaresl Il98sl) that the apparent luminosities 
of many galactic nuclei are several orders of magnitude 
smaller than is predicted by Bondi-Hoyle accretion onto 
the black holes that are known to reside there, given 
plausible lower limits on the central gas density in the 
cooling flow. The Advection Dominated Accretion Flow 
(ADAF) model of accretion onto black holes has been exten- 
sively explored a s a way of explaining this surprising result 
Jlchimarul Il977l : iNaravan fc Yilll994T) . However, the basic 
physical premise of the ADAF model, that coulomb scat- 
tering is the dominant process for establishing cquipartition 
between electrons and ions, is probably false <IBinnevl 2003). 
iBlandford fc Begelmanl (l999) propose an alternative to the 
ADAF model as an explanation of the low luminosities of 
galactic nuclei. In their ADIOS model most of the accretion 
energy that is released as plasma falls into the black hole is 
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used to drive a wind from the surface of the accretion disk. 
Most of the gas that falls onto the outer edge of the accre- 
tion disk is carried by this wind away from the black hole, 
with the result that the hole's accretion rate is much smaller 
than the disk's accretion rate. 

Our simulations assume that the ADIOS model cor- 
rectly predicts that much of the energy released by accretion 
onto the black hole is channelled into a sub-relativistic bipo- 
lar flow. We emphasize that this flow is probably distinct 
from the highly relativistic flow that generates jets of syn- 
chrotron radiation in many well observed systems, such as 
M87. Key differences between the bipolar flow on which we 
focus and the underlying synchrotron jets include (i) the for- 
mer is subrelativistic since material is blown off the accretion 
disk at the local Kepler speed, while superluminal motion 
and one-sidedness show the latter to have a Lorentz factor of 
several; (ii) the former will comprise an ordinary H/He gas, 
while the latter is quite possibly an e± plasma; (iii) the for- 
mer is a direct and inevitable corollary of accretion while the 
latter probably draws it s energy from the black ho le's spin 
via vacuum breakdown jBlandford fc ZnaiekHlirnl) . and it 
may switch on and off in an erratic way. 

The distinction between the bipolar flow from the disk 
and the one that generates synchrotron jets is important be- 
cause the latter flow probably dominates radio maps at all 
radii through its superior ability to generate ultrarelativis- 
tic electrons. However, in many observed systems both flows 
are likely to be present simultaneously, so we should imagine 
the outflow to consist of a series of concentric cylinders (or 
cones in the case of non-negligible jet opening angle), each 
cylinder moving parallel to the axis at a speed that increases 
from ~ 100 kms" 1 on the outside to ~ c at the centre. Insta- 
bilities powered by the shear within this system will steadily 
transfer momentum and energy outwards, and cause ambi- 
ent plasma to be entrained at the edge. A very basic point, 
but one that is often overlooked, is that the mechanical lu- 
minosity of the entire jet is orders of magnitude greater than 
the synchrotron luminosity. The ratio of these luminosities 
will be especially large if there is no ultrarelativistic jet at 
the core of the subrelativistic outflow. 

Outflows with speeds ~ 0.1c and mechanical luminosi- 
ties of order the Eddington luminosity are di rectly observed 
in sp ectra of accreting relativistic objects iPounds et al. I 
2003). While these outflows are suggestive, they are not di- 
rectly relevant to the systems of interest here because they 
are observed in objects with compara ble photon luminositie s 
and they may be radiatively driven jKing fc Poundsll2003l) . 
It is likely that in very dense environments much of the ac- 
cretion energy is degraded into photons, while in systems 
with very hot and therefore rarefied central gas, the accre- 
tion energy emerges in largely mechanical form. 



4 SIMULATION TECHNIQUE 

Even with the best current soft- and hard-ware it is im- 
practicable to simulate flows from the sub-parsec scales on 
which jets form, out to beyond r coo i ~ 100 kpc. In practice 
one must start at some smallest resolved scale r m j n ~ 1 kpc, 
using a model of the jets on that scale that derives from a 
mixture of physical intuition and a critical examination of 
observed jets. From observations and modelling of jet evo- 



lution at r ~ r m i n , one hopes to infer the structure that the 
jets must have on much smaller scales if they are to have 
the inferred structure of scales r m i n and above. 



4.1 Generating jets 

On a given scale, a jet is characterized by the rates m, P and 
E at which it injects mass, momentum and energy into the 
larger-scale plasma. If the characteristic speed of flow within 
the jet, wjet, is subrelativistic, these three numbers are clearly 
related by P = mvjet and E — m(u + |w? et ), where u is the 
specific internal energy of the jet material. So long as the 
jet is highly supersonic, u <C ^vf ct . Entrainment causes rh 
to rise and Vj e t to decline as one proceeds down the jet, such 
that P remains approximately constant. Entrainment leads 
to an increase in u relative to the kinetic contribution to E, 
while radiative cooling causes u to decrease. 

Our simulations employ the three-dimensional 
adap tive-mesh hydrocode ENZO fervan fc Norman] 
Il997t) with a Piecewise Parabolic M ethod (PPM) Riemann 
solver dColella fc Woodwardl Il984h . The refinements are 
dynamically generated based on the density gradient down 
to 6 levels of refinement, with an effective resolution of 
1024 3 cells. The computational box is 635 kpc on a side, 
and as a result the central ~ 35 kpc is covered by cells 
measuring 620ft" 1 pc. Periodic boundary conditions are 
enforced. 

We simulate the action of a pair of sub-grid jets by 
adding mass, i-momentum and energy to cells that lie 
within the yz plane at x = ±e, where e < 1 kpc. During 
a timestep of length 8t, the injected mass rhSt is distributed 
over cells that lie in a thin disk according to the window 
function 

w{r) oc\ e Jct for r < Hot, (1) 

[ otherwise. 

Here r = */ y 2 + z 2 is the distance of the cell from the x- 
axis and fjet = e/0.3. The injected momentum and energy 
are distributed in the same manner. We present results for 
two values of 7>t , namely 2 and 3 kpc to illustrate the effect 
of changing this parameter. The cross-sectional area of the 
jet at the highest resolution is resolved into ~ 35 (78) cells 
for the rjct = 2 kpc (3 kpc) case. 

The first increments in the dynamical variables after a 
jet is switched on increase the mean velocity of material in a 
disk cell only very slightly because the injected momentum 
has to be shared with a substantial quantity of stationary 
gas. Consequently, little of the injected energy is tied up in 
the kinetic energy of the cell. Hence, nearly all the injected 
energy is used to heat the gas that was originally in the cell, 
and this gas expands in response. The resulting decline in 
the density of gas in disk cells causes subsequent momentum 
increments to yield ever larger bulk velocities for disk cells, 
and gradually a significant fraction of the injected energy 
goes into bulk kinetic energy. 

From this discussion it will be seen that the speed Vj e t 
that characterizes the ratio of the increments in momen- 
tum and energy, controls the relative importance of heating, 
which tends to expand the nuclear gas isotropically, and bulk 
motion, which generates bipolarity within the cooling flow. 
In this paper we restrict ourselves to simulations in which 
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Figure 1. The period of small radial oscillations as a function of 
radius in the initial model. 



Vjet = 10 000 km s~ and m = 2Mq yr~ , which corresponds 
to a total power P — 6 x 10 43 erg s _1 for two jets. The jets re- 
main powered for 100 Myr so that they inject 2 x 10 59 erg in 
total. We shall explore the effects of changing the jet power 
in a later paper. The low value of Vj et adopted here contrasts 
with the much higher let ve l ocitie s (~ 100 00 km s -1 ) stud- 
ied in iRevnolds etafl feoolEool and lBasson fc Alexander! 
J2003I) . 

Our approach to jet simulation avoids the imposition 
of boundary conditions at points internal to the simulation. 
It involves adding source terms for each cell to generate the 
appropriate m, P, and E. Given these, the code computes 
increments in the internal energy u. Neither the tempera- 
ture nor the density at the base of the jet is hard-wired; 
they evolve as natural consequences of m, P, and E. Our 
approach guarantees a gradual dynamical response to the 
jets being either turned on or off. 

The power of our jets is smaller than the X-ray lumi- 
nosities of the clust er that we model (Lx = 2.5 x 10 44 erg s _1 
JPavid et all200(fl 1. so jets of this power could significantly 
modify the cooling flow only if they were active for much of 
the time. It seems likely that jet activity is intermittent, and 
our goal in this paper is to understand the impact of a sin- 
gle outburst. The relevant dynamical timescale of the cluster 
gas turns out to be quite long (see below), so the dynamical 
impact of an outburst is evident more than a gigayear after 
the jets have switched off. By contrast, in the absence of a 
further outburst, the core will proceed to a cooling catastro- 
phe within ~ 300 Myr fe.g. lKaiser fc Binnevl i200StO . Hence 
if we are to understand the impact of a single outburst, we 
must suppress radiative cooling. We reserve to a later pa- 
per simulations that include radiative cooling and repeated 
outbursts. 



4.2 The cluster model 

The background that is disturbed by the jets is based on 
Chand ra observations of the Hydra cluster by iDavid et alJ 
(2000). The gravitational potential is that of the NFW 
model that David et al. fitted to their data. From this poten- 
tial we calculated electron densities for an isothermal distri- 
bution of plasma, and chose the temperature T = 3.7 x 10 7 K 
that provides the best fit to the electron density that David 
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Figure 2. Temperature and x-velocity as a function of time for 
points lying on the jet axis at x = 2, 12 and 45 kpc. 
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Figure 3. Electron density and the x component of velocity as 
functions of perpendicular distance from the jet axis sampled at 
x = 12 kpc. Figs. HI and 151 show two-dimensional slices through 
the density field at each of the times plotted here. 



et al. inferred. The equations of motion were integrated for 
100 Myr before the jets were switched on to ensure that the 
system was not far from numerical equilibrium when the jets 
were ignited. The jets fired for 100 Myr before turning off, 
and the subsequent evolution was followed for 2.6 Gyr. 

The jets deform the initially spherical surfaces of con- 
stant specific entropy. After the jets have turned off, the 
surfaces of constant specific entropy are expected to wiggle 
in and out with a characteristic period that increases from 
small radii to large. The period of these oscillations are of 
order r = 2ir/ui, where the Brunt- Vaisala frequency lo is 
given by 



2 d$dln<7 

u> — 



dr dr 



(2) 



where $ is the gravitational potential and a = Pp' 1 is the 
entropy index. In Fig. we plot r as a function of radius 
in our initial model of the Hydra cluster. The period rises 
from ~ 0.2 Gyr at the origin to 1 Gyr at 150 kpc. Hence the 
period is everywhere longer than the likely lifetime of a jet. 
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Figure 4. Density in the plane z = at 6 to 117 Myr after the jets turn on, in the case of the 2kpc jets. The intensity scaling and length 
scale is the same for each image, but the last two panes show only the positive x-hemisphere. 



5 RESULTS 

Fig- El describes the process of jet formation once mass, mo- 
mentum and energy start being injected into the disk by 
plotting the temperature and z-velocity on the jet axis at 2, 
12 and 45 kpc from the cluster centre. The upper full curves 
show that the temperature at x = 2 kpc rises steeply to a 
sharp peak as injection starts. The temperature then fluc- 
tuates before settling to a plateau that lies above 10 9 K in 
the case of the narrower jet, and slightly below 3 x 10 s K in 



the wider case. The dotted curves in the temperature panels 
show that after a delay this process is repeated at slightly 
lower temperatures at x = 12 kpc. The full curves for the 
x- velocity in Fig.|2]show that at x = 2 kpc the velocity rises 
slightly less rapidly than the temperature before falling and 
levelling off at a 4500 kms -1 plateau for the narrow jet, and 
1800 km s" 1 for the wider one. Since the speed of sound in 
the ambient medium is 920kms~ 1 , the corresponding Mach 
numbers with respect to the ambient medium are M = 5 
and M = 2. The jets' internal Mach numbers are much less 
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Figure 5. Same images as in Fig. |1] but for the case of the 3kpc jets. 



than unity. For both jets the peak velocities are larger at 
x — 12 than at 2 kpc. In light of the lower temperatures at 
12 kpc, this suggests that between these two points the jet 
has narrowed and accelerated as in a nozzle. Thus these are 
pressure-confined jets such as are thought to exist in low- 
power (FR I) radio galaxies, rather than the ballistic jets of 
FR II sources. 

The dot-dashed curves in Fig. [5] show similar trends at 
x — 45 kpc except that (i) the impact of the jet is delayed 
by ~ 45 Myr, and (ii) the temperature attains much smaller 
peak values. The lower temperatures again suggests adia- 



batic expansion. The pressure at a given radius should be 
approximately time-independent since it is determined by 
the ambient cooling flow. Hence, in the absence of entrain- 
ment, the Bernoulli function of each fluid element would be 
conserved. Actually the decrease in T is not associated with 
a corresponding increase in velocity, so entrainment proba- 
bly is important. Direct measurements of mass flux at vari- 
ous radii confirm this inference. At all three locations both 
temperature and velocity drop precipitately soon after in- 
jection ceases. Interestingly, the delay between the velocity 
dropping from location to location is smaller than the corre- 
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sponding delay that accompanied jet ignition. Thus it takes 
only 15 — 20 Myr for the jet channel to fill up out to 45 kpc, 
compared to the ~ 45 Myr required to cut the channel that 
far. 

The fact that the plateau velocities at x — 2 kpc in 
Fig- El are substantially smaller than i>j e t reflects the impor- 
tance of entrainment at the base of the jet. Near the plane 
x = 0, gas is sucked in towards the x axis, flows through 
the injection disks, and is then blasted out along the axis. 
This flow diminishes the final jet velocity in two ways: (i) it 
diminishes the pressure just behind the injection disk, which 
tends to throttle the flow; (ii) it causes the momentum in- 
jected on the disk to be shared by a larger quantity of gas. 
The flow induced through the back of the injection disk is 
larger in the case of the wider jet, because the area of its disk 
is 2.25 times that of the narrower disk. In fact, the plateau 
velocity of the wider jet is smaller than that of the narrow jet 
by a very similar factor, 2.5. The detail of this entrainment 
process is highly artificial, but the general physical principle 
is sound: gas will be entrained by a jet that has propagated 
out from the AGN to the radius at which we have placed 
our injection disks. By measuring the flux of gas at 2 kpc 
we find that the 2 kpc jet entrains ~ 3M0 yr -1 at its base, 
while the 3 kpc jet entrains ~ IOM0 yr -1 . 

Fig. |3 shows for five times electron density and x- 
velocity along a cross-sections through the jets at x = 
12 kpc. (The density fields in a two dimensional plane at 
these times are depicted in Figs. 0] and |SJ) The upper pan- 
els in Fig. |21 show that by 26 Myr, injection has drastically 
lowered the density on the axis, by a factor ~ 25 in the case 
rjet = 2 kpc and by a factor ~ 6 in the case of the wider 
jet. The situation at 18 Myr is more complex: in the nar- 
rower jet the density is already low out to ~ 4 kpc of the 
jet axis, but it is slightly raised in the region beyond 4 kpc 
by the jet-driven flow of material away from the centre; at 
this early time the wider jet has yet to reach x = 12 kpc so 
the density lies above its original value throughout the cross 
section and, at 300 km s -1 , the velocity is still subsonic. At 
18 Myr the narrower jet is highly supersonic (with respect 
to the ambient medium; v ~ 6000 kms -1 ) and shows a sig- 
nificant backfiow in the region y = 2 — 5 kpc, with speeds 
up to 970 kms -1 . The wider jet develops a backfiow later 
at x = 12 kpc, and in both cases the backfiow region moves 
away from x = 12 kpc at later times. It is interesting that 
both the density and velocity plots show that the jets be- 
come narrower with time. This fact reflects the movement 
downstream of the backfiow region and suggests that the 
artificial shear viscosity in ENZO is rather small. The ve- 
locity profiles for t — 55 and 76 Myr show quite wide wings 
either side of the narrowing core. Although the velocities as- 
sociated with these wings are modest, the corresponding net 
momentum is a significant fraction of the jet's momentum 
because the areas and densities associated with the wings 
are large. Indeed, the directly measured mass flux in the re- 
gion 2 < 2/ kpc < 12 of either jet reveals a further mass 
loading of 8 - 1OM yr" 1 . 

Figs.2]and|S]show the density field in a two dimensional 
plane (z = 0) at the same six times for the 2 and 3 kpc jets, 
respectively. At the first time shown, i.e. 6 Myr after the jets 
were turned on, two lunes of low density can be seen, one on 
each side of the injection disk. These lunes rapidly swell par- 
allel to the x axis and develop internal structure that makes 
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Figure 6. Density as shown in the t = 117 Myr panel of Fig. IS1 
but with an altered intensity scaling which enhances the equiden- 
sity surfaces. 



them resemble mushrooms, those for the wider jet being 
shorter and plumper. By 55 Myr more complex structures 
have developed that differ significantly between the jets: the 
narrower jet has a butterfly-shaped cavity about 15 kpc long 
that has its head about 50 kpc from the origin and a wide 
tail around a long straight section of jet that reaches 30 kpc 
back to the origin. By contrast, the cavity generated by the 
wider jet extends only to x ~ 40 kpc along the axis and is 
shaped rather like a red pepper. A shorter, wider jet leads 
back from this structure to the origin. The structures seen 
at t = 76 Myr show similar differences between the two jets. 
At 117 Myr, 17 Myr after the jets turned off, the channels 
along which the jets ran are rapidly filling near the origin 
and the outer structures are becoming more turbulent. The 
heads of the cavities continue to plough outwards. About 
100 Myr after the jets have turned off (not shown) the out- 
ward crashing structures undergo a fundamental change: a 
stream of colder than ambient gas surges up from behind 
the cavity, largely filling it and bursting out of its leading 
edge. In this way the cavities turn into overdensities. 

In the 76 Myr snapshot of the 2 kpc jet we see that the 
flow is disturbed in the region 25 < x/ kpc < 35, where the 
jet no longer follows a linear path. This is a highly turbulent 
region with strong internal shocks that disrupts the jet stem 
and almost separates it from the established cavity. There 
is no analogous feature in the 3 kpc jet, but do we find that 
such features play an important role in simulations with 
faster jets. 

Fig. |S] which is a centered view of the cluster core at 
t — 117 Myr, shows that the equidensity surfaces remain 
markedly aspherical some time after the stem has retreated 
from a given radius. They soon after become roughly ellip- 
tical and distort under radial oscillations that are analogous 
to the g-modes of a star. Fig. [7] quantifies the relaxation of 
the equal-entropy surfaces by plotting as a function of scaled 
time t/r, where r is the local Brunt Vaisala time from Fig.0 
the quadrupole moment of the entropy distribution 
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Figure 7. Quadrupolc moment of the density distribution as a 
function of t/t s v with t < 2.6 Gyr at r = 25, 50, 110 and 200 kpc. 
Curves for successive radii have been shifted horizontally by 0.03 




-wo -100 o 100 



Figure 8. Maps of specific entropy in the plane z = at late 
times. 

Q(r,t)= ddsmO , (3) 

7 2 <r(r,*J 

where 9 is colatitude with respect to the jet axis and W(r, t) 
is the average value of the entropy index a on a sphere of 
radius r. At r = 25 and 50 kpc a large excursion in Q is 
evident as the jet first impacts, followed by an oscillation 
pattern that looks like an interference pattern between the 
first few harmonics with fundamental period comparable to 
r(r). At r — 110 and 200 kpc very few cycles are seen even 
though the baseline in time extends to 2.6 Gyr. 

Fig- |U shows the entropy field for the 3 kpc jet at three 
late times. The panel for 380 Myr shows the flow after the 
cavities have transformed into over-densities. In the panel 
for 710 Myr we see that around x — 70 kpc overdense mate- 
rial has fallen back and indented the otherwise fairly smooth 
isoentropy surfaces at that radius. In the bottom panel we 
see the same process happening further out at 1 Gyr. This 



outward-moving wave of material falling back strongly ex- 
cites the system's g-modes. 

ENZO follows the density of several dyes, each one of 
which moves with the main fluid according to the usual con- 
tinuity equation. This facility enables us to explore the ex- 
tent to which an outburst moves material radially. The full 
curves in Fig. [0] show the distribution at six times of a dye 
that was initially distributed like the gas density interior to 
r = 5 kpc, and had zero density outside this 'inner core' vol- 
ume. Data for the narrower jet are on the left. The dashed 
curves show the distribution of a second 'outer core' dye, 
which was initially distributed like the gas density in the ra- 
dial range 5 < r/kpc < 77 and had zero density elsewhere. 
The upper middle panel for the narrower jet in Fig.[5]shows 
that 47 Myr after the jet turns on, a significant amount of 
inner-core dye has moved out to the range 20 < r/kpc < 40; 
the next panel shows that when the jet turns off, the inner 
core material extends out to 100 kpc. The lower panels show 
that after the jet has died, inner core material is carried out 
to r > 200 kpc before falling back slightly in the last part of 
the simulation. Similar effects are evident in the right panels 
for the 3 kpc jet, the main difference being that this wider 
jet ejects essentially all rather than most of the inner core 
material. 

The dashed curves in Fig.[5]show that the jet has a less 
dramatic impact on the outer-core material. Nevertheless, in 
the left panels we see that outer-core material is pushed out 
to beyond r — 100 kpc, mostly after the jet has switched off. 
Moreover, by the end of the simulation with rj e t = 2 kpc, 
sufficient outer-core material has flowed in to the inner-core 
region for the densities of the two fluids to be essentially 
equal at the centre, with the density of the outer-core ma- 
terial declining much more slowly with radius than that of 
the inner-core material. The equality of the central densi- 
ties is masked in Fig. [5] by the difference in the adopted 
mass scales. The ability of the wider jet to displace most 
of the inner core gas while affecting the outer core gas only 
about as much as the narrower jet does, is consistent with 
the conclusion we reached from Fig. |3] that the main differ- 
ence between the two jets lies in the different quantities of 
plasma that they entrain near the injection discs. 

Fig. 1101 shows the evolution of the metallicity of the 
IGM. This was followed by distributing a dye in the ini- 
tial configuration with a density pd that is proportional to 
a power of the main plasma density, with the power cho- 
sen such that the 'metallicity' Z — pd/ P declines by a fac- 
tor 10 between r = and r = 500 kpc. The upper set of 
four curves, which shows the spherically averaged metallic- 
ity density for the simulation with the narrower jet, shows 
that the expulsion of inner-core material weakens the central 
metallicity gradient slightly, but does not eliminate it. The 
wider jet causes the central metallicity to decline slightly 
more: to 0.857 rather than 0.888 of its original value. Natu- 
rally changes in metallicity are largest along the line of the 
jets. The lower set of curves in Fig. 1101 shows as a function 
of distance x down the jet axis the metallicity an observer 
would measure along the jet if the cluster were oriented such 
that the jets lie in the plane of the sky. The original metal- 
licity gradient is smaller than that obtained by spherical 
averaging. Also at x ~ 20 the metallicity is clearly seen to 
be enhanced by uplift of material from the core. 

The top row of Fig. Illl shows the divergence of the ve- 
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Figure 9. The distribution in log radius of material that started at r < 5kpc (full curves) or at 5 < r/kpc < 77 (dashed curves). The 
vertical scale is in units of 5 X 10 8 Mq for the full curves, and 5 X 10 10 Mq for the dashed curves. Results for the narrower (rj ot = 2kpc) 
jet are shown at the left with results for the jet with rj ct = 3kpc on the right. 
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Figure 10. Evolution of the metallicity of the cluster gas in the 
simulation with the narrower (2kpc) jet. Upper curves: spheri- 
cally averaged metallicity; lower curves: metallicity derived from 
projected data at points along the projected jet axis. 



locity field in the plane z = at 14 and 55Myr after ig- 
nition of the 2kpc jet. At the first time shown we can see 
the injection disks and, nearly connected to them, a network 
of bright regions that roughly mark the boundaries of the 
mushrooms we encountered in the density field. Further out 
a circle of brightness delineates the bow shock as it advances 
into the undisturbed IGM. At 55 Myr the divergence field 
is more complex, on account of the emergence of internal 
shocks within the jets. A pair of wedges mark the interface 
between ejected material and the disturbed IGM and a little 
further out a bow shock appears as two sectors of a circle. 
The original quasi-spherical shock is too weak to see in this 
plot, but the unsharp- masked X-ray image at bottom right 
reveals it at about 10 o'clock and 2 o'clock. At 3 o'clock the 
jet's bow shock has just penetrated the spherical locus of the 
original shock. The distance between the jet head and the 
original shock has not increased between 14 and 55 Myr, so 
the head is moving supersonically - from Fig. 2]we find that 
the jet head moves 37 kpc in the 37 Myr from t = 18 Myr, 



so its average speed is a factor 1.08 times the sound speed 
at 3.7 x 10 7 K. 

The middle row in Fig. 1111 shows the X-ray emissivity 
projected along the z axis at the times of the top panels. 
The shocks that are clearly evident in the top panels are in- 
visible in the X-ray data. The bottom panels show unsharp- 
masked versions of the mid dle panels. In these im ages the 
shocks are evident. Recently iFabian et al] <l2003al) have un- 
sharp masked a very long exposure of the Perseus cluster 
and detected numerous weak arc-like regions of enhanced 
density gradient that are similar to those seen in Fig. 1111 

The right-hand panels in the middle and lower rows of 
Fig. llll show X-ray cavities that are very si milar to those ob- 
serve d i n Hydra, Perseus a nd other clvjgte^sjjBoIrnn ger et_al J 
[ 19931: iMcNamara et all l2000t IFabian et al. I l200oF 
Heinz et al. I 120021 iBlanton. Sarazin fc McNamaral I2003D . 
The cavities in our simulations are slightly edge brightened, 
but rather less so than are the observed cavities because the 
material at the edge is not cooler than ambient. In part this 
will be because we start with an isothermal temperature 
profile rather than one that drops near the centre. It may 
also be the case that in a realistic cluster environment the 
jets plough into dense material dredged up by a previous 
jet. This would produce a bow shock comprised of denser 
material, which would shine more brightly in X-rays. It is 
therefore important to extend these sim ulations to include 
multiple jet eruptions dSoker et al.ll2002T) . 

Fig. H2l shows the overall impact that the jet has on the 
system by comparing the radial run of mass at four times for 
the narrower jet (the corresponding plots for the wider jet 
are similar). The upper panel shows the radius r(M) within 
which mass M lies. This increases slightly at the centre. The 
lower panel shows this increase in more detail by showing 
the ratio of the radius r t that contains M at time t to the 
radius ro that originally contained the same mass. Near the 
centre, r t is ~ 20 percent bigger than ro so long as the jet 
is powered. Further out a wave of increase in rt can be seen 
propagating out while the jet is on. By 350 Myr after the 
jet is switched off, r t is larger than ro by < 5 percent at all 
radii. 

We have seen that the jet transports cool material from 





Figure 11. Top row: the divergence of the velocity field in the plane z = at 14 and 55 Myr after jet ignition shows the location of the 
bow shocks. Middle row: the X-ray emissivity projected along the z axis at 14 and 55 Myr show little or no trace of the shocks. Bottom 
row: the shocks can be clearly seen in unsharp-masked versions of the X-ray images of the middle row. In the 55 Myr case, the central 
region has been blocked out to enhance the weak bow shock. All images are for the rj et = 2 kpc case. 



the centre and deposits it at large radii. Hence an increase in 
r(M) at fixed M does not necessarily mean that a physical 
shell of material moves outwards; the shell may in fact move 
inwards but its decrease in r be more than compensated by 
the decrease in its mass coordinate M. 

Fie. ll3l shows the distribution of binding energy per unit 
mass as a function of mass at four times scaled such that 
the area under the curves is total binding energy in units of 
10 58 erg. Thus the injection of 2 x 10 59 erg by the jets has 
the effect of shifting the mean level of the curves vertically 
by only 0.1 in the units of the vertical axis. The lower panel 
shows the difference in the binding energy per unit mass at 



each enclosed mass between the given time and the moment 
before jet ignition. Heating by the jet reduces binding energy 
and thus lowers the curves. Around M = 40 x 10 10 Mq we 
see a pronounced dip in the curve for t — 47 Myr. There is 
a corresponding dip in the curve for t — 101 Myr around 
M = 130 x 10 10 M Q . This increase in the radius of the dip 
reflects outward propagation of energy. The decline on the 
extreme right in the curve for t = 447 Myr implies that by 
this time the extra energy has migrated to the highest mass 
shells included in the figure, which lie near r = 200 kpc, 
appro ximately the cooling radius in this cluster jDavid et alJ 
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Figure 12. M(r) at four times for the narrower jet. 
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Figure 14. Entropy index versus mass at four times for the nar- 
rower jet. 
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Figure 13. |dE/dlnM| at four times (units 10 58 erg per 
1O 1O M0) for the narrower jet. 



Consider the effect of uplift of cold gas on a physical 
shell of matter that lies between the radii from which the 
jet takes and then dumps cold gas. If the shell is not di- 
rectly disturbed by the jet, its energy will be constant, but 
the value of its mass coordinate M will decrease. Since in 
the original configuration specific binding energy |d£/dM| 
decreases with increasing M, the decrease in M will have the 
effect of decreasing |di?/dA/| at given M. Hence in Fig. 1131 
uplift has the same effect as a distributed source of heat that 
acts between the radii between which the jet transports ma- 
terial. 

The areas under the curves of the lower panel should 
equal the energy injected up to each time and thus provide 
valuable checks on the accuracy of the simulation. The area 
under the full curve should be 9.4 and is 9.45, while the area 
under the other two curves should be 20 and is 22.6 for t = 
101 Myr and 11.4 for t — 446 Myr. Extremely similar values 
are obtained for the jet with rj ct = 3 kpc. In view of the 
smallness of the energy injected relative to the total energy, 
the agreement between the expected and measured values at 
the first two times is satisfactory. The shape of the dashed 
curve in Fig. 1131 suggests that the discrepancy at the third 
time arises because the integration has not been carried out 
to sufficiently large enclosed masses - the outermost mass 



shell included lies at r ~ 200 kpc. Thus this figure indicates 
that the energy injected by the jet ends near the cooling 
radius. A simple calculation explains why this is so. We have 
seen that the jet blows a cavity that travels at close to the 
sound speed and remains an underdensity for about twice 
as long as the jet fires. Hence the radius at which it becomes 
an overdensity and deposits its energy is ~ 2c s ij c t ~ 180 kpc 
for our choice tj ct = 100 Myr. If heating is to make good 
radiative losses, the heating action of the jet needs to be 
concentrated at rather smaller radii. Thus our assumed jet 
lifetime is likely to be too large by a factor of a few. Since 
the size and energy content of the cavities we generate are 
realistic, any decrease in tj ct should be compensated by an 
increase in jet power. 

Fig. 1141 shows the variation with enclosed mass of the 
entropy index a = Pj 1 p 1 in the case of the narrower jet. A 
prompt increase in the entropy density at M <C 10 x 1O 1O M0 
is evident. Subsequently, the entropy density decreases in 
this zone, but remains above its initial value. However, 
around M = 25 x 10 10 M Q the entropy index at t = 446 Myr 
drops slightly below its initial value. This depression is asso- 
ciated with an enhancement in binding energy that is clearly 
visible towards the left-hand edge of Fig. 1131 From Fig. ll2l wc 
see that the corresponding radial coordinate is r ~ 30 kpc, 
which Fig. 1111 shows to be the extent down the x axis of 
a grid boundary. Consequently, there has to be a suspicion 
that this feature is a numerical artifact. 



6 COMPARISON WITH OTHER WORK 

Churazov^^dJ 




iQuilis et al] fcOOll) . 
iBriiggen et al. I (120021) and 

injected energy only, whereas 



Briiggcn & Kaiser 
Briig gen fc Kaiser! 
we have injected both energy and momen tum. Moreover, 
with the exception of IQuilis et al.l i200lf) these authors 
simulate only a part of the clust er that varies i n scale 
from the 30 x 10 2 kpc volume of |Bruggen et al. I (120021) 
to a hemisphere of radius 1 M pc of |BruggeriJ^£a iseil 
<200ll). Only the sim ulations of IQuilis et al] feOOll) and 



Briiggen et al~l <2002ft are fully three-dimensional, and only 



Briiggen fc Kais er (2002) had an adaptive grid: unfortu- 



nately, this two-dimensional simulation had slab rather 
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than rotational symmetry. In all these simulations energy 
was injected at a fixed point that is significantly removed 
from the cluster centre without any physical motivation for 
its location. 

The presence of the bow shock (see Fig. II It and of the 
strong turbulence in the ejected material underlines the dif- 
ference between our model and models in which 'bubbles' 
of plasma buoyantly rise up in a quasi-stationary manner 
IChurazov et alJl200ll; iQuilis et alJl200ll : IrBriiggen fc Kaiserl 
120011. 120021: iBrtiggen et al. II2002I) . Outward motion of the 
point at which the jet's energy is thermalized is clearly im- 
portant for the dynamics of cavities in our simulations, and 
missing from most early discussions of bubble dynamics. 

Two groups have simulated the impact of jets on cool- 
ing fl ows (iRevnolds et alj|200ll. 120021 : iBasson fc Alexander! 
2003). Both groups employed a version of the ZEUS code 
on a single spherical grid that excluded a sphere around 
the origin. In the case of Reynolds et al., the radius of this 
sphere was 1 percent of that (500 kpc) of the outer boundary 
sphere, while in the simulations of Basson & Alexander the 
ratio of these radii was 0.5 percent and the outer bounding 
radius was 2800 kpc. Our computational domain is 630 kpc 
on a side and includes the origin. In these simulation the 
PPM Riemann solver was employed, while ZEUS uses arti- 
ficial viscosity to simulate shocks. 

While Reynolds et al. forced axisymmetry, the simula- 
tions of Basson & Alexander were fully three-dimensiional. 
In both simulations the jets were imposed as conical inflows 
through inner boundary. The jets were faster and more pow- 
erful than ours: they emerged from the inner sphere at a 
speed of ~ 100 000 km s -1 and a density 100 times lower 
than the ambient density. In the simulations of Reynolds 
et al. the net jet power was 9.3 x 10 45 ergs _1 . We esti- 
mate the power of the jets of Basson & Alexander to be 
2.5 x 10 46 ergs _1 . Hence, in both simulations the jets are 
more powerful than ours by at least a factor 150 to 400. 
The rate of mass injection was just 1.5 to ~ 4 times ours. 
Thus the jets were significantly less heavily loaded with mo- 
mentum than ours, yet there was enough momentum for the 
head, at which energy was randomized, to move out through 
the cluster, rather than remaining stuck at an arbitrarily 
chosen point, as in earlier work. 

An important difference between the Reynolds et al. 
simulations and ours is that in their simulations, the back- 
flows extended all the way to the origin for the lifetime of 
the jet, whereas in our simulations the backflow region moves 
away from the origin at some distance from the head of the 
jet. This difference is a natural consequence of the different 
momentum loadings of the two jets since a higher ratio of 
momentum to energy increases the speed of advance of the 
head relative to the speed of the backflow. 

A second notable difference between the simulations is 
in the temporal and spatial scales of the outbursts. In the 
Reynolds et al. simulations, the cavities created by the jets 
have still not been overtaken by uplifted cold material at 
15 times the lifetime (50Myr) of the jets, by which time 
the cavities are ~ 1 Mpc from the centre. Thus in these 
simulations the AGN has a major impact to the edge of the 
cluster. A similar situation is evident in Fig. 1 of Basson & 
Alexander. The observations suggest that the AGN's impact 
should lie primarily within the cooling radius and in future 
work we plan to concentrate on the part of parameter space 



in which this is so. The high-power outbursts simulated by 
Reynolds et al. and Basson & Alexander clearly do not have 

promising parameters from this point of view. 

Al l simulation s exce pt t hose of |Bjajj£gen_fc_i£ajse]J 

ll2002ft . lQuilis et al] fcOOlft and lBasson fc Alexander! (I2003D 
neglect cooling, as we do. 



7 DISCUSSION 

It is now clear that AGN have a significant impact on cool- 
ing flows. The largest question that remains unresolved is 
whether the best observed systems are in an approximate 
steady state, or are cycling between configurations that have 
significantly different central density profiles and tempera- 
tures. Another important issue is whether the observed cavi- 
ties are driven outwards by momentum-laden jets, or simply 
rise buoyantly. We have argued that the low luminosities of 
black holes embedded in X-ray emitting gas, coupled with 
the im plausibility of the ADAF model, support the conjec- 
ture of iBlandford fc Begelmanl <ll999T) that a subrelativistic 
bipolar flow is the primary output channel of these black 
holes. This conclusion amounts to a prima facie case for 
momentum driving. Our simulations explore the regime of 
strong momentum driving, and provide a counterpoint to 
simulations in which energy alone is released at some point 
in the flow. 

We have modelled the impact on a cooling flow similar 
to that in the Hydra cluster, of an outburst by an AGN 
that lasts 100 Myr and has mechanical power 6 x 10 43 erg s _1 
delivered by opposing jets. The jet power is smaller than 
the X-ray luminosity of the Hydra cluster, so to achieve a 
quasi-steady state a series of more powerful outbursts would 
have to follow one another. We have sought to gain physical 
insight into the overall response of the cooling flow to energy 
input by jets by suppressing both subsequent outbursts and 
radiative cooling. 

The jets are created within the simulation by an algo- 
rithm that avoids the imposition of any boundary conditions 
near the cluster centre. The speed and early development 
of the jets depends on the width of the base of the jets 
because wider jets entrain more ambient material at their 
base. The gross features of jet action are independent of 
the assumed jet width, however. Cavities of very hot plasma 
grow around the heads of the jets, and these move super- 
sonically outwards. Ahead of the cavities run two cap-like 
bow shocks, one over each jet head, which eventually over- 
take the weakening spherical shock front that is generated 
as the jets ignite. The cavities, which advance supersoni- 
cally, are evident in simulated raw X-ray images, but the 
bow shocks can only be detected in the X-ray data when 
they are unsharp-masked. 

A turbulent vortex that contains a significant quantity 
of entrained and uplifted material trails each cavity. After 
the jet has turned off the cavity weakens and slows relative to 
the denser vortex. About 100 Myr after the jet turned off, the 
vortex fills and overtakes the cavity, which is thus replaced 
by a region of outward-moving, cool, overdense gas. Eventu- 
ally this over-den sity falls back inwards. Many studies (e.g. 
lAlexanderl j2002l) and references therein) have pointed out 
that an AGN outburst will lift significant quantities of low- 
entropy material out of the cluster core. However, they have 
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discussed this process in the context of shock-acceleration of 
material ahead of a rising cavity. It is not clear that shock- 
acceleration causes the dense following vortex that we see; 
this structure is probably a part of the global circulation 
pattern that the jets excite, in which gas flows in perpen- 
dicular to the jet axis and out parallel to the axis. 

We have used dyes to assess the importance of turbulent 
mixing. Much or most of the material that was originally in 
an inner core of radius 5kpc is ejected to radii as large as 
200 kpc. Material that originally lay between 5 and 77 kpc 
moves in part inward to replace the ejected inner-core ma- 
terial, and in part outwards. Most stays put, with the re- 
sult that an original metallicity gradient in the IGM is only 
slightly weakened. 

There is a small overall inflation of the IGM: while the 
jets are still firing, the radius at a given value of the mass 
coordinate M increases by ~ 10 percent within a few kilo- 
parsecs of the centre. The long-term increase in radius of 
material within ~ 100 kpc is only of order a percent. The 
entropy index of the innermost ~ 1O 11 M0 increases by up 
to a few percent. Most of the energy injected by the jets 
ends up around the cooling radius or beyond because the 
distance that the cavities move is of order 2c a tjet. Heating 
could more readily offset cooling if a smaller value of tj e t and 
a larger value of the jet power were used. 

We now return to the questions posed in Section 3: 

(i) The outflow from an accreting black hole must carry 
momentum as well as energy. Most previous studies have ne- 
glected the momentum flux because they have been directed 
at the effects of ultrarelativistic jets, that probably consist of 
a very light electron-positron plasma. We have stressed the 
likelihood that any such jet lies at the core of a massive, sub- 
relativistic bipolar flow that comes from the accretion torus 
rather than the black hole's ergosphere. Our assumption of a 
non-negligible momentum flux enables our simulations to be 
self-consistent in a way that the energy-only simulations are 
not: in our simulations the point at which the jet's ordered 
kinetic energy is randomized is self-consistently determ ined, 
whereas in the simulations of lBriiggen fc Raised (l2002t> and 
others the point at which energy is injected is arbitrarily 
chosen. Moreover, cavities that are pushed out into the am- 
bient medium are more stable than ones that merely rise 
buoyantly and we tentatively infer from the morphology of 
the Virgo A radio source that the data favour stable cavi- 
ties. Specifically, Fig. 1151 shows that a featur e in th e radio- 
continuum map of Virgo A in lOwen et al.l ll200Cf) closely 
resembles a feature in one of our simulations. This resem- 
blance cannot count as hard evidence for momentum driving 
because the quantities plotted in the two cases are very dif- 
ferent. But our impression is that no feature in the energy- 
only simulations so closely resembles the observed feature 
because the simulated features are too prone to fragmenta- 
tion. 

(ii) The jet's energy offsets the effects of radiative cool- 
ing in several ways: 1) gas at the very centre is pushed up 
high, where it mixes in with higher-entropy material and 
is heated by dissipation of the system's strongly excited g- 
modes; 2) gas slightly further out is first shock heated and 
then adiabatically compressed as it sinks inward to take up 
the volume vacated by ejected gas; 3) gas at all radii is 
heated by being turbulently mixed with the extremely hot 




Figure 1 5. Left: a piec e of the radio-continuum map of Virgo 
A by iOwen et alj|2000l) . Right: a piece of the map of specific 
entropy in the plane z = from one of our simulations. 

shocked jet material; 4) the radio map of lOwen et all l l200Ch 
strongly suggests that in Virgo shocked jet material is dis- 
tributed throughout the inner ~ 40 kpc. If the jet consisted 
primarily of ultrarelativistic material, it would not follow 
that the thermal plasma had gained significant heat from 
the jet because the coupling between ultrarelativistic parti- 
cles and thermal plasma is weak. If we accept the likelihood 
that the observed ultrarelativistic jet is sheathed in a sub- 
relativistic wind, it does follow from the radio map that tur- 
bulent mixing will have heated the thermal plasma within 
~ 40 kpc of the cluster centre. 

(iii) The rate of jet heating can be estimated as the en- 
ergy within the observed cavities, divided by cavity lifetime 
~ 2tjet. Our conclusion that cavities survive for about twice 
the time it takes to inflate them is consistent with the ob- 
servation of two pairs of cavities in clusters such as Perseus, 
one pair being radio bright and the other a pair of ghost 
cavities. Since cavities move at the speed of sound, we can 
determine their lifetimes from how far out we find them. 
For Perseus we thus infer ij e t = 15 Myr, and a heating rate 
P — 1.4 x 10 44 ergs -1 per 7kpc-radius cavity. This is close 
to but less than the X-ray luminosity from within the cool- 
ing radius. So are cooling flows in steady states? We think 
not. First, it is unlikely that the black hole's output can be 
closely matched to the radiative cooling rate because the 
characteristic timescales of the bla ck hole and the co oling 
flow are so discr epant. Second, as iMotl et al. I feOOSt) and 
iKav et ail i2003l) have stressed, clusters experience substan- 
tial infall of gas, and a significant quantity almost certainly 
reaches the centre without being shock heated to the virial 
temp erature. It is likely that the gas observed in mole cular 
form ( Donahue et al. l200dlEdgd200ltlEdge et all2002h and 
as Hq-emitting filam ents iLyndaTHiTOnCbwie^tHn. Ill983l : 
IConselice et aljfeoOlF) is ju st such gas rather than gas that 
has co oled in the flow: as ISparks. Macchetto fc Golombekl 
( 1989) have stressed, dust in the filaments would be unlikely 
to have perfectly normal Galactic extinction properties if the 
gas had condensed from the hot phase. Moreover, the fila- 
ments would not have their highly disturbed morphology if 
they had condensed from the flow. X-ray data from Chandra 
show that the col d gas is in thermal co ntact with the X-ray 
emitting plasma jFabian et alj|2003rj) . and it is probably 
being evaporated by the latter. If so, this activity will be a 
significant drain on the energy resources of the cluster gas, 
making it likely that, notwithstanding current AGN heat- 
ing, clusters such as Perseus are drifting towards a cooling 
catastrophe. 
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(iv) After an outburst, hydrostatic equilibrium is ap- 
proached after several Brunt- Vaisala times. From FigJ^we 
see that this timescale is very long, being > 1 Gyr at r — 
200 kpc and still larger further out. Consequently cooling- 
flow clusters should show significant deviations from hydro- 
static equilibrium at most radii, and several different AGN 
outbursts will have an impact at the current time. 

(v) Raw X-ray images are expected to show cavities but 
not shocks, even though the cavities move supersonically. If 
X-ray data of exceptionally high signal-to-noise are avail- 
able, unsharp-masking will reveal the shocks. Cavities live 
for about twice as long as the jet that inflates them, so there 
should be 2-4 cavities visible at any given time. We do not 
simulate the production of synchrotron-emitting electrons, 
but in view of the short lifetimes of these particles (~ 40 Myr 
for Lorentz factors 7 ~ 10 4 ), one would not expect cavities 
to be bright synchrotron sources for long after their parent 
jet died. Hence our simulations suggest that a significant 
fraction of cavities should be radio-invisible 'ghost' cavities. 

(vi) Although Fig. 1141 shows that the specific entropy of 
matter at the centre is increased by the AGN, Fig. 1131 sug- 
gests that most of the energy injected by an outburst ends 
up near or beyond the cooling radius. This result follows 
from the facts 1) that much of the jet's energy goes into 
inflating a cavity which moves mildly supersonically, and 2) 
that we have assumed that the jet fires for 100 Myr and find 
that the cavity survives for a further 100 Myr after the jet 
expires. Travelling at ~ 1000 km s for 200 Myr, the cavity 
moves out to ~ 200 kpc before depositing its energy. 

If one accepts that averaged over gigayear times, heating 
by AGN offsets radiative cooling, then the heating ought 
to be concentrated well within t he cooling radius, which is 
~ 20 kpc in the Hydra cluster JPavid et al.ll200dl . More- 
over, [Owei^^ii] (120001) find that in Virgo the synchrotron 
emission has a sharp outer boundary at ~ 60 percent of the 
cooling radius. This observation suggest s that in Virgo no 
cavity has reached the cooling radius, and lSoker et alJ <l2002l) 
remark that observed radio jets do not propoagate to large 
distances in cooling-flow clusters. We tentatively conclude 
that in our simulations the cavities move to unrealistically 
large radii. This deficiency would be resolved by halving the 
lifetime of the outburst. Since the energy in our cavities is 
about right, this halving of the lifetime would have to be 
offset by a doubling of the jet power. 

(vii) A crude picture of the impact of jets on the ICM 
is that the jets entrain a substantial body of gas from the 
very centre, heat it and transport it to rather large radii, 
where it mixes in with ambient material. Material that orig- 
inally overlay the ejected gas then sinks inwards to take its 
place. In this picture the curves of metallicity versus radius 
in the left panel of Fig. llOl would simply stretch to the left at 
constant metallicity because the gas from the centre would 
have a negligible effect on the metallicity of the much larger- 
quantity of gas with which it mixed at large radii. In this 
case the decline in the central metallicity gradient would 
be very small. Fig. 1101 indeed shows that, although an out- 
burst flattens the metallicity gradient, the effect is confined 
to r < 20 kpc and is small. Ongoing stellar evolution would 
have a countervailing tendency to increase the metallicity 
gradient. Thus there does not seem to be a conflict between 
heating by jets and the existence of metallicity gradients 
near cluster centres. 



(viii) Since Fig. shows that the dynamical time of the 
ICM exceeds 200 Myr at all radii, the large-scale structure 
of the ICM will be unaffected by fluctuations in the AGN's 
power on short timescales. However fluctuations on shorter 
timescales may be important. Of particular interest is the 
timescale on which the jet channel closes up after the jet 
turns off, since a jet that is quiescent for this time or longer 
can probably not reactivate the cavities that it was formerly 
inflating. Hence new cavities will be formed when the jet 
comes back online. Our simulations suggest that the chan- 
nel fills up more quickly than they are cut, but shorter 
timescales are conceivable, especially if the ICM is rotating 
and the jet does not fire along the axis of rotation. 

Although great strides have been made in our under- 
standing of cooling flows in the last few years, many key 
questions remain open. The cooling flow phenomenon is an 
extremely important one both because of its connection with 
galaxy formation, and because of the insights it gives to the 
dynamics of AGN. Simulations will have an important role 
to play in understanding cooling flows because they are in- 
herently dynamical systems. Although a certain amount of 
progress can be made with axisymmetric simulations, the 
chaotic and turbulent nature of cooling flows makes essen- 
tial three-dimensional simulations such as those we have pre- 
sented here. Limitations of the simulations that arise from 
only limited spatial resolution are a worry, particularly in 
the three-dimensional case. For this reason we have cho- 
sen to focus in this paper on just two simulations, which 
are highly artificial in that they exclude radiative cooling. 
These simulations have clarified a number of numerical is- 
sues and provided a foundation for further work. We are 
currently experimenting with shorter-lived, more powerful 
jets and with simulations that include cooling and rotation 
of the cluster gas. Animations of the simulations described 
here and a number of other simulations can be viewed at 
www . clust erheat ing.org. 
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